I'm not really sure that this is needed. You could just get the border with a P_1 finite element space and then multiply and shift to get a "multi-component" border.

I'm not really sure that this is needed. You could just get the border with a P_1 finite element space and then multiply and shift to get a "multi-component" border.

You are right. I ignored the shortcut.

But for [P2, P2, P1], this seems to be more complicated, Could I do it following your suggestion in the case?

It won't be as easy for [P2, P2, P1], but your code won't work either for such a fespace, I guess.

I try it, from the result of plot, it looks like ok when I run the following code, and the number of DOFs on the border seems ok.

include "ExtractDofsonBorder.idp"

macro fe() [P2, P2, P1] //

mesh Th = square(64, 64);
int[int] labs = 1 : 3;
int[int] doflabs;
fespace Vh(Th, fe);
ExtractDofsonBorder(labs, Vh, doflabs, 1, 3);
Vh [uh, vh, ph] = [0, 0, 0];
for (int j = 0; j < doflabs.n; ++j) {
    uh[][doflabs[j]] = 1;
}
plot(uh, cmm="uh", wait=1);
plot(vh, cmm="vh", wait=1);
plot(ph, cmm="ph", wait=1);